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Abstract 

Homeostatic control of cell volume and intracellular electrolyte content 
is a fundamental problem in physiology and is central to the functioning 
of epithelial systems. These physiological processes are modeled using 
pump-leak models, a system of differential algebraic equations that de- 
scribes the balance of ions and water flowing across the cell membrane. 
Despite their widespread use, very little is known about their mathemat- 
ical properties. Here, we establish analytical results on the existence and 
stability of steady states for a general class of pump-leak models. We 
treat two cases. When the ion channel currents have a linear current- 
voltage relationship, we show that there is at most one steady state, and 
that the steady state is globally asymptotically stable. If there are no 
steady states, the cell volume tends to infinity with time. When minimal 
assumptions are placed on the properties of ion channel currents, we show 
that there is an asymptotically stable steady state so long as the pump 
current is not too large. The key analytical tool is a free energy relation 
satisfied by a general class of pump-leak models, which can be used as a 
Lyapunov function to study stability. 

Key Words: Cell Volume Control, Electrolyte Balance, Free Energy, Lya- 
punov Function, Differential Algebraic System. 



1 Introduction 

Cells contain a large number of organic molecules that do not leak out through 
the cell membrane. The presence of organic molecules and their attendant coun- 
terions results in excess intracellular osmotic pressure. The plasma membrane 
is not mechanically strong enough to withstand significant differences in os- 
motic pressure, and thus the cell will tend to swell and burst. Plant cells and 
bacteria have a mechanically rigid cell wall to guard against this tendency. An- 
imal cells maintain their cell volume with ionic pumps and ionic channels that 
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together regulate the ionic composition of the cytosol ( Hoffmann et al. , 20091 



Evansl . l2009t iBoron and Boulpaepl 120081) 



This "pump-leak" me c hanism is typically modeled in t he following fashion 
( Keener and Snevdl . Il998t iHoppensteadt and Peskinl . 120021 ). Consider a cell of 
volume v and let be the intracellular and extracellular ionic concentrations 
respectively. We only consider the ions Na + , K + and CI - . Let the cell be in 
a large and well-stirred extracellular bath so that [-] e can be assumed constant. 
We have the following balance equations for the three ionic species. 
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Here, </> is the membrane potential, g. are the ion channel conductances for 
each species of ion, a is the strength of the pump current, F is the Faraday 
constant, and RT is the ideal gas constant times absolute temperature. The 
pump current for Na + and K + has a ratio of 3 : 2 reflec ting the stoichiometry of 
the Na-K ATPase (in (jHoppensteadt and Peskinl 120021 ). this ratio is set to 1 : 1 
for simplicity). The above balance laws are supplemented by the following: 
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Equation (jl.ldp is the electroneutrality condition, where A is the total amount 
of organic molecules in the cell and z the average charge on one organic molecule. 
We have assumed that there are no organic molecules outside the cell and that 
they do not pass through the membrane. Equation (|l.lej) says that water flows 
into or out of the cell according to the osmotic pressure difference across the 
membrane. Here, f is the membrane permeability to water flow. System (jLTj 
forms a system of differential algebraic equations. We would like to see under 
what condition the above system possesses a stable steady state, representing a 

cell with a stable cell volume and ionic compositiom 

Models of t his type were first introduced in ( Tosteson and Hoffman . 19601 : 
Tostesorl 1964 ) and have since been extended and modified by several authors to 
study cell volume contro l (jjakobsson . Il98d: iLew et aTlll99ll : lHernandez and Cristinal . 
1998c lArmstrong) . 120031) . Pump-leak models are widely used in epithelial phys- 



iology. Epithelial cells need to maintain their cell volume and ionic composi- 
tion in the face of widely varying extracellular ionic and osmotic conditions. 
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There is a large body of mathematical modeling work in the context of renal 
physiology (see (jWeinsteinl . Il994 120031 ) for review). Mat hematical modeling 



studies of other systems us i ng pump - leak m odels include (jLarsen et al 
iFischbarg and Dieckei [200l lYi et all 120031 ). 



2002 



Despite their widespread use and fundamental physiological importance, 
there seem to be very few an alytical results regarding t he behavior of pump- leak 
models. As for system (jl.ip . (jKeener and Snevd . 1998t) shows the following. As- 
suming z < — 1, there is a unique steady state with a finite positive cell volume 
if and only if: 



[Na+] c cxp(-3ai7( ffN ai?T/F)) + [K+] e exp(2aF / (g^RT / F)) 



[Na + ] 
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(1.2) 



(1.3) 



is satisfied, system (jl.lj) possesses a steady state for sufficiently small a > 0. 

To the best of the author's kno wledge, there are n o analytical results on the 
stability of these steady states. In ( Weinsteinl . Il997l ) the author studies an ep- 
ithelial model of greater complexity than (jl.lj) . The author obtains an algebraic 
expression for the linearized matrix around steady state and numerically stud- 
ies its eigenvalues for physiological parameter values. The computation of this 
linearization is complicated by the presence of an algebraic constraint (the elec- 
troneutrality condition). Some authors have considered simpler non-electrolyte 
models of (epithelial) cell volume con trol and solute t r ansport. Analytical re- 
sults for such models can be found in (|Weinsteinl . Il992l : iHernandezl . 120031 . 120071 : 
" 2010h . 



Benson et al 



The goal of this paper is to establish analytical results on the existence and 
stability of steady states for a large class of pump-leak models that includes 
(jl.lj) and other representative models as special cases. In Section [21 we intro- 
duce the general class of pump-leak models that we shall treat in this paper. 
We consider iV-species of ions subject to the electroneutrality constraint. Cell 
volume is controlled by the transmembrane osmotic pressure difference. The 
key observation of this Section is that the system of equations satisfies a free 
energy identity. This identity and its variations will be the main tool in study- 
ing the stability of steady states. The presence of a free energetic structure in 
pump-leak models leads us to natural structure conditions to be imposed on 
current- voltage relationships for ionic channel currents. 

In Section [3l we study the case when the ionic channel current (or pas- 
sive ionic flux) has a linear current-voltage relationship, the pump currents are 
constant and water flow is linearly proportional to the transmembrane osmotic 
pressure difference. System (jl.lj) is an example of such a system. We first es- 
tablish the necessary and sufficient condition under which the system possesses 
a unique steady state. This condition, when applied to system (jl.lj) . reduces to 
(jl.2j) (in fact, our conclusion is slightly stronger; we shall see that the restriction 
z < — 1 is not needed in (|1.2p ). We then prove that this steady state is globally 
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asymptotically stable. If a steady state does not exist, the cell volume v tends to 
infinity as time t — > oo for any initial condition. The main tool in proving these 
statements is a modified version of the free energy introduced in Section [5] This 
modified free energy G satisfies dG/dt = — J, J > 0, thus defining a Lyapunov 
function. Asymptotic stability of steady states follows by an examination of 
G. To prove the global statements, we also make use of the fact that J is a 
Lyapunov function. This is a consequence of the fact that, in a suitable set of 
variables, the system is a gradient flow of the convex function G with respect to 
a suitable metric. We thus have a clear dichotomy; if the system has a steady 
state, it is globally asymptotically stable and if not, the cell bursts. In Section 
13.31 we discuss a simple epithelial model in which the cell is in contact with a 
mucosal and serosal bath. When the current voltage relationships for the ionic 
channels are all linear, the same stability results hold for this simple epithelial 
model. 

In many modeling studies using the pump-leak model, the current-voltage 
relation for the ionic channel current is not linear. The popular Goldman cur- 
rent voltage relation is one such example. The goal of Section @] is to establish 
a result on the existence and stability of steady states with minimal assump- 
tions on the current-voltage relation. Indeed, all we assume here are properties 
required on thermodynamic grounds. We first discuss solvability of the differ- 
ential algebraic system. Solvability is not entirely trivial given the algebraic 
constraint of electroneutrality. We then show that for a general class of pump- 
leak models, there is a steady state for sufficiently small positive pump rates so 
long as a generalization of condition (|1.3[) is satisfied. We then show that this 
steady state is asymptotically stable. Our main tool here is again the modified 
or relative free energy. The main difficulty in establishing stability is that the 
current voltage relation cannot in general be written as a function of the chem- 
ical potential jump only (this difficulty is not present when the current voltage 
relation is linear). We shall see that this effect is small when the pump rate is 
sufficiently small, which allows us to establish asymptotic stability of the steady 
state. 

2 Model Formulation and the Free Energy Iden- 
tity 

Consider N species of ion and let Cfc, k = 1, • • • , N be the intracellular concen- 
trations of the fc-th species of ion. We let c e k be the extracellular concentrations 
of these ions, which are assumed to be positive and constant independent of 
time. Let v be the volume of the cell. The balance equation for the ions can be 
written as follows: 

j t (vc k ) = -Mc,c c )-p k ( ( f>,c,c ), k=l,---,N. (2.1) 

Here, we have written the transmembrane flux as the sum of the passive flux 
jk and the active flux pk- The active flux, typically generated by ionic pumps, 
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requires energy expenditure whereas the passive flux, carried by ionic chan- 
nels and transporters, does not. The flux functions jk and pk depend on the 
transmembrane potential <f> (intracellular potential minus the extracellular po- 
tential) as well as the vector of intracellular and extracellular concentrations 
c = (ci, ■ ■ • , cm) t and c c = (c°, • • • , c e N ) T ', where - T denotes the transpose. We 
assume that jk and pk are C 1 functions of their arguments. Since the c c arc 
assumed constant, the c° only appear as parameters of the above differential 
equation. The dependence on c e will thus often not be shown explicitly. 

The functional form of the active flux function pk can be arbitrary, but for 
jk, we must impose structure conditions so that it represents a passive flux. 
Commonly used examples of jk are: 

j\ = G k (iff In f ^ J + FzkA = G k ^ k , Gk > 0, (2.2) 



■GHK _ p 

Jk — r k 



Fz k <f> C k ex P (tH^ 
^ exp(^)-l 



exp \ RT 



where i?T is the ideal gas constant times absolute temperature, F is the Faraday 
constant and Zk is the valence of the /c-th species of ion (e.g., 1 for Na + , —1 for 
Cl~ and so on). We assume that there is at least one ionic species for which 
Zk 7^ 0. In the above, fik is the chemical potential of the A:-th intracellular ion 
measured with respect to the extracellular bath: 

Hk =RT In (^] +Fz k 4>. (2.4) 



Expression j\ is linear in \ik- If we multiply this by Fzk to change units from 
ion flux into electric current, we obtain a linear current voltage relationship for 
this ionic current. This was used in Expression jj? HK can be derived by 

assuming a constant electric field across the ion channel, and is known as the 
Goldman-Hodgkin-Katz current formula. An important feature of both j\ and 
jj? HK is that they are increasing functions of fik for fixed <f> and that it is when 
\ik is 0. In the general case, we require jk to satisfy a somewhat weaker version 
of these properties, which we shall discuss later in relation to Proposition ^. II 

Observe that Ck can be expressed in terms of fik and <j> through (|2.4p . We 
shall often find it useful to view jk as a function of cj> and fi = (//i, • • • ,^n) T 
instead of </> and c. We shall write this as jk{4>, A 1 ) m a slight abuse of notation. 
We note that the dependence of jk on /x; , I ^ k expresses the possibility that 
the flow of the k-th ion may be driven by the chemical potential gradient of the 
Z-th ion. This is the case with many ionic transporters in which the flow of one 
species of ion is coupled t o another. Indeed, many models of io nic transporter 
currents have this feature ( Weinstein . 1983 ; Strieter et al. . 1990l ). 
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Equations (|2.ip are supplemented by: 
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(2.5) 
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j w (c,c e ,v). 



(2.6) 



In (|2.5j) , A > is the total amount of organic molecules inside the cell, and z is 
the average valence of intracellular organic molecules. Equation (j2.5p states that 
both the intracellular and extracellular concentrations satisfy the electroneu- 
trality constraint. Electroneutrality of the extracellular space, together with 
the requirement that z k ^ for at least one k, requires that there must be at 
least two ionic species. In (|2.6p . j w is the passive transmembrane water flux. 
An example of j w is: 



where £ > is the hydraulic conductivity of water through the membrane. This 
simple prescription is what is used in (jl.lj) and many other models of cell volume 
and electrolyte control. In this example, water flow is proportional to 7r w , the 
osmotic pressure difference across the membrane, whose expression is given by 
the van't Hoff law. We shall assume that j w is a C 1 function only of 7r w . The 
important property of (|2.7[) is that j w = when 7r w = and that it is increasing 
in 7r w . This property will be discussed further in relation to Proposition 12. II 

We seek solutions to the differential algebraic system (|2.1j) . (|2.5|) and (|2.(il) 
given initial values for c = (ci, ••■ ,cn) t ,v and <j) that satisfy the algebraic 
constraint (|2.5|) . We require that Cft > 0, k = 1, • • • ,N and v > for all time. 
The membrane potential <j> evolves so that the electroneutrality constraint (|2.5p 
is satisfied at each instant. Multiplying (|2.ip by Fz k and summing in k, we 
have: 



where we used (|2.5p to conclude that /, the total transmembrane electric current, 
must be 0. This gives us an algebraic equation for <p that must be solved at 
each instant. The solvability of this equation is a necessary condition for the 
initial value problem to be solvable. For certain specific functional forms of jk 
and pk, like those of (|2.2p and (|2.3p with p k constant, the solvability of (|2.8[) is 
immediate. We shall discuss the general case in Sectional 

One way to avoid the above difficulty arising from the algebraic constraint 
is to replace (|2.5p with the following relation for intracellular concentrations: 




(2.7) 




(2.8) 




(2.9) 



G 



Mathematics of Pump-Leak Models 



Mori 



where C m is the total membrane capacitance. This states that the total charge 
inside the cell is equal to the charge stored on the membrane capacitor. We note 
that (|2.9p has its own difficulties as a biophysical model. If the concentrations 
are defined as the total amount of intracellular ions divided by volume, (|2.9[) 
is indeed correct. If we define Ck to be the ionic concentration away from the 
surface charge layer (on the order of Debye lengt h sa lnm in width), we must 
introduce surface ionic densities as was done in dMori and PeskinL l2009h . As 



we shall see shortly, the left hand side of (|2.9p is often negligibly small and 
the electroneutrality condition (|2.5|) can thus be seen as a perturbative limit of 
condition (|2.9|) . We shall discuss this point further after we make our system 
dimensionless. 

Scale ionic concentration, volume and membrane potential as follows and 
introduce the primed dimensionless variables: 

RT 

Cfe = c c' k , c c k = c c k , v = v v', <j) = —0', (2.10) 

where cq and vq are the typical concentrations and volumes respectively. Equa- 
tions (|2"TTj) . and (|2l?|) become: 

d -^ = -f k { (j) '^')- P ' k ^,c'), ^=ln(^+W, (2.11a) 

JL A' * A 

dv' N ( N A'\ 

% = -fMi < = £ car - £ 4 + £ ), (2.iic) 
fc=i \fc=i / 

where /x' = (fi[ , ■ • • , n' N ) T and c' = (c[ , • ■ ■ , c' N ) T . Time t and the flux functions 
jkiPk and j w are suitably rescaled to yield their respective primed variables. We 
note that it is possible to further reduce the number of constants, for example, 
by taking cq to be the total extracellular concentration. We shall not pursue 
this here, since it leads to some difficulty in understanding the physical meaning 
of each term in the resulting dimensionless system. 
Equation (|2.9[) yields: 




e = C ™ RT ' F (2.12) 

FcqVq 

where e is a dimensionless parameter expressing the ratio between the amount 
of ions contributing to the surface charge and the absolute amount of charge 
in the cytosolic bulk. This quantity is typically very small (about 10 -7 ) and 
we thus expect that it is an excellent approximation to let the left hand side of 
(|2.12p be and adopt condition (|2.5p (or its dimensionless version (|2.11bp ). if the 
membrane potential does not vary too rapidly. Most, if not all modeling studies 
of cellular electrolyte and water balance use the electroneutrality constraint 
(|2.5p or (|2.11bp and we shall treat this case only. 
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We shall henceforth deal almost exclusively with the dimensionless system. 
To avoid cluttered notation, we remove the primes from the dimensionless vari- 
ables. 

An important property of the above system is that it possesses a natural 
energy function. 



Proposition 2.1. Let v andck,k 
following equality holds: 



1, • • • , N satisfy system (|2.11[) . Then, the 
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(2.14) 



Identity (|2.13p is not entirely new. In (jSauer , 1973; lFromterl .11974; Wei nsteinl . 

1983h the authors argue on thermodynamic grounds that free energy dissipa- 



tion and input for an epithelial system should be given by the right hand side 
of (|2.13|) . The new observation here is that, by appropriately defining a free 
energy function (that is to say, the left hand side of (|2.13[) '). this thermodynamic 
property can be turned into a mathematical statement about system (|2.11[) . We 
also point out that a similar identity valid for a system of partial differential 
equat ions describing electrodiffusion and osmosis was proved in ( Mori et all 
l20lfl . In subsequent sections we will use this as a tool to study stability of 
steady states. 

Proof of Provosition [Ql View a as a function of c k and ca = A/v. Note that: 

d 



dt 



(vc A ) 
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~dt 
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since A is constant, 
molecules as: 



Note that: 



Define the chemical potential of intracellular organic 

da 

dc A 



HA = In c A + zq 



zqb. 



(2.16) 



^ V z k <p- 



(2.17) 

Multiply (|2.1[) by fi/~, multiply (|2.15[) by ji A and take the sum. The left hand 
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side yields: 

N d d 
5>* 5 («*)+A*a 5 («u) 
fc=i 

/c— i x 

d . . / <9cr <9cr \ cfo <i / s \ 



d dv 
= dt {va) -^ 

where we used (|2.5[) and ()2.14|) in the third equality. We thus have: 
d dv N 

— (va) - 7T W — =-^2 +Pk)- (2.19) 



k=l 



Equation (|2~T3| thus follows from (|2~6| . □ 

The function a should be interpreted as the free energy per unit volume of 
intracellular electrolyte solution. A key fact that was used in the above proof is 
the identity: 

(v-^ da da \ 

g^ + ^J' (2 - 20) 

This relation, connect ing the free energy with osmotic pressure, is well-known 
in physical chemistry (jDoi and Seel . ll996h . 

When pk = 0, fc = 1, • ■ ■ ,N in (|2.13|) . there are no active currents and we 
have: 

dG N 

-TT = ~^2^kjk ~ TTwJw (2.21) 
k=l 

Given the interpretation of G as the total free energy of the system, the second 
law of thermodynamics requires that G be decreasing in time. The negativity of 
(|2.2ip when /i ^ and 7r w ^ is equivalent to the statement that the following 
conditions be satisfied: 

JV 

M3k(<l>, M) > for all and fi ^ 0, (2.22) 



k=l 



7r w i w (vr w )>0ifvr w ^0, (2.23) 

Condition (|2.22[) , together with continuity of jk and j w with respect to its 
arguments immediately implies that: 

j*(&A* = 0) = 0, k = !,-■■ ,N, j w (7r w = 0)=0, (2.24) 
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where the conditions on jf. is to be satisfied for all (f>. Taking the derivative of 
the above expression for jk with respect to <fi, we see that: 

||(&Ai = O)=0, k = l,--- ,N. (2.25) 

We shall find this expression useful later on. 

Let us require that the derivative of j w with respect to 7r w be non-zero at 
7r w = 0. This non-degeneracy condition, together with (|2.23l) . leads to: 

f^K - 0) > 0. (2.26) 
oir w 

Let j = (ji, ■ ■ ■ , jm) t and let dj/dfi be the Jacobian matrix with respect to fx 
for fixed <f>. That is to say, the kl entry of the N x N matrix dj/dfi is given by 
djk/dm. The non-degeneracy condition for j is that dj/d/i be non-singular at 
fi = 0. We require the following condition: 

d\ 

■77— (6, u = 0) is symmetric positive definite for all 6. (2.27) 

Of! 

The symmetry of the Jacobian matrix does not follow from the non-degeneracy 
condition and condition (|2.22p . In fact, these two conditions imply only that: 



0) is symmetric positive definite for all 4>. (2.28) 




The symmetry of dj/dn is required by the Onsager reciprocity princip le ( Onsager 



1931t iKatzir-Katchalskv and Curranl . Il965t iKielstrup an d Bedeaux, 2008). We 



note that (12.271) i s the same as the condition introduced in (|Sauerl . ll973tlFromter 



lci pie [\um 
uxl.l2008l). 



1974 IWeinsteinl . 1 19831 ). It is easy to see that both $TZ§ and (gUJ satisfy flUg) , 
(|2~24]l and ([2^27)1 and that (gj]) satisfies ([2~23]) . (|2T24)) and (|2~26l) . 

Note that (|2.23p together with (|2.11cp implies that the intracellular and 
extracellular osmotic pressures must be equal at steady state. We are thus 
assuming that the membrane cannot generate any mechanical force to balance a 
difference in osmotic pressure. If the cell membrane (or its attendant structures) 
can generate some elastic force, it would make it easier for the cell to maintain 
its volume. 

Our starting point in deriving the above structure conditions for jk and j w 
was the requirement that the right hand side of (|2.21[) be negative. If jk is 
allowed to depend on ir w and j w on /x, a more general structure condition can 
be formulated. Although it should not be difficult to extend the results to follow 
to this more general case, we will not pursue this here to keep the presentation 
reasonably simple. 
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3 Results when the Flux Functions are Linear 
in the Chemical Potential Jump 

Before dealing with the general case in Section 0] we treat the simpler case when 
jk and j w are linear in fj, and 7r w respectively and pk are constants in (12. lip . 
In this case, we obtain a more or less complete picture of the behavior of our 
system. System (|2.11[) becomes: 

4(«c) = -Lfi-p, (3.1a) 

n zA N 

= VzfeC fe + — =V z k c%, (3.1b) 

i < y Z / 

k=l k=l 
Tt = (3 - 1C) 

where L = dj/d/j, is an TV x TV matrix, which by (|2.27[) , is symmetric positive 
definite, and p = (pi, • • • ,Pn) t is the vector of the active fluxes. The hydraulic 
permeability £ = dj w /dTr w is positive by (|2.26[) . The extracellular ionic concen- 
trations c| > 0, k = 1, • • • , TV and the amount of impermeable organic solute 
A are assumed positive as discussed in the previous Section. We seek solutions 
(c, v) G x M + where R + denotes the set of positive real numbers. 

Here and in the sequel, we shall often find it useful to refer to the pair 
(c,i>) as well as the triple (c,v,<fi) G x R + x E. We shall often view the 
pair and the triple as being members of and x R respectively and 

write (c,v) e and (c,v,<T>) £ R^ +1 x M C M. N+2 . The more "correct" 

notation may be to write (c T ,v) T = (ci,-- - ,cn,v) t g (c t ,v,(/)) t — 

(ci, ■ • • ,cn,V, (f>) T G x M given that c is a column vector. We will not 

adopt this unnecessarily ugly notation. Similar comments apply to the pair 
(a, v) and the triple (a, i;, <ft) where a = vc. 

For system (|3.ip . we may compute 4> explicitly in terms of c by solving the 
(dimensionless version of) (|2.8I) : 

<t>= / — r \ ' T= (Ti)-" ,7jv) ,7fc=ln — , (3.2) 

^z, i^z; K jv \ c fc/ 

where z = (zi, • • • ,zjv) T - Substituting the above expression for </> into p. lap , 
we obtain: 

d i \ 7 1 , ^ ? t (Lz) k (Lz)i 

— (vc ) = -L(7 + q), L kl =L kl - , 

di (z,£z) R iv (3.3) 

q = (si, • • • ,?jv) t , q = i _1 p, 

where Lki,Lki are the fcZ entries of the TV x TV matrix L, L and (iz)^ is the 
fc-th component of the vector Lz G K w . Note that L -1 exists given that L is 
positive definite. Replacing p. lap with (|3.3p . we now have a system of ordinary 
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differential equations (ODEs) for c and v only. Applying the standard existence 
and uniqueness theorem for ODEs, we conclude that a unique solution always 
exist for sufficiently short time so long as the solution remains in (c, v) € R^ +1 . 

Even though our ODE system is for N + l variables c = (ci, • ■ ■ , cn) t and v, 
our dynamical system is only N dimensional, since the dynamics is constrained 
by the electroneutrality condition (|3.1b[) . The initial value problem for p.l[) can 
thus only be solved if the initial values satisfy ()3.1bp . In the proof of Proposition 
13.41 we will find it useful to make a change of variables to remove this constraint. 



3.1 Existence of Steady States and Asymptotic Stability 

Our first observation is the following. 
Proposition 3.1. Consider the function: 

N 

f((ft) = ]T c% (exp(-q k - z k (ft) - 1) , (3.4) 
fc=i 

where q k ,k = 1, • • • ,N was defined in (|3.3p . The function f(<fi),4> € K has 
a unique minimizer (ft = (ft m i n . System (|3.ip has a unique steady state if this 
minimum value is negative: 

imin(q, c c , z) = f(<ft min ) < 0. (3.5) 

Otherwise, the system does not have any steady states. 

The above condition can be interpreted as follows. At steady state, it is easily 
seen that the concentrations c k must be equal to c£ = c c k exp(— q k — Zk<ft*) where 
(ft* is the value of (ft at steady state (see (pTT)) ). We need f((ft*) = J2k=i( c t~ c k) < 
since there must be "osmotic room" for the impermeable solutes. This is only 
possible if the minimum of /(</»), (ft G K is negative. We also point out that the 
above condition depends only on q = L~ 1 p, c c and z and does not depend on z 
or A. 

Proof of Provosition \3. 1\ Set the right hand side of (|3.1a[) and (|3.1c[) to zero. 
We have: 

q = £~ 1 p = -M, 7r w = 0. (3.6) 
Solving for c& in the first expression we have: 

c k = c% cxp(-q fc - Zk<j>), k = 1 • • • ,N. (3.7) 

Substitute this into ir w — and (|3.1bj) . We have: 

(3.8) 
(3.9) 





A 
f — 

V 


= 0, 




4f + 

d(ft 


zA 

V 


N 

= ^2 z kC c k cxp(-q k - z k (ft) 
fc=i 


zA 
+ — =0 

V 
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where f(<j>) is given by (|3.4p . We must find solutions <fi and v > to the above 
system. Note that: 

d 2 f ^ df 

"77^ = E z l c k cx P(-% - z k4>) > 0, lim -L = ±00. (3.10) 

^ kl 0±oo dip 

The second property comes from the fact that there are ions with negative and 
positive valences among the N species of ions and that c| > 0. We can thus 
solve (|3.9|) for <p uniquely in terms of v. Let this function be <f> = ip(v). We have: 



dip ( d 2 f\ 1 zA 



(3.11) 



dv 

Consider the left hand side of (|3.8j) and substitute <p — ip(v) into this expression: 

R(v) ee f(<p(v)) + ^ = 0. (3.12) 

Our problem of finding steady states is reduced to the question of whether the 
above equation in v has a positive solution. We have: 

dR = tfd^_ A = _(flY 1 M!_ A <_A <0 (3 13) 

dv dip dv v 2 \ dcf> 2 J v 3 v 2 ~ v 2 

where we used (|3.9[) and (|3.11[) in the second equality and (|3 . 10[) in the first 
inequality. Therefore, R(v) is monotone decreasing. Note that: 

R(e) = R(l)- £ f^j dv > R(X) + J (Jp) dv = i?(l) + A(e- 1 -l). (3.14) 

Therefore, R(v) — > 00 as v tends to from above. Thus, (|3.12l) has a unique 
positive solution if: 

lim R(v) < 0, (3.15) 

and otherwise, there is no solution. Let ip^ = lim^oo ip(v). Note that this 
limit exists since, by (|3.11[) , ip(v) is monotone if z and constant if z = 0. 
Taking the limit v — > 00 on both sides of (|3.9[) . we see that cpoo is the unique 
solution to df /dip = as an equation for <p. Given (|3.12j) . condition (|3.15l) can 
be written as /(y>oo) < 0. The statement follows by taking m ; n = ipoo. □ 



Condition (|3.5p applied to yields condition (|1.2p . Since condition (|3.5p 
is valid regardless of the v alue of z, we may lift the restriction z < — 1 found in 
([Keener and Snevdl . Il998h . 

Fix q = L x p and c c so that ()3.5p is satisfied. Since condition p.5p does 
not depend on z or A, a unique steady state (c,v,<fi) = (c*,v*,<fi*) exists for 
any value of z and A > 0. We may thus view (c*, v* , <p*) as a function of A and 
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Qa = zA, defined for A > and Qa £ M. We can compute the dependence of 
v* on Qa and A as follows. 



dv* 
~dA 

dv* 




d 2 f 



>0, 



(3.16) 
(3.17) 



From a biophysical standpoint, (|3.16[) is reasonable since more impermeable so- 
lute leads to greater osmotic pressure. Note that (|3.17l) says that the v* increases 
if the absolute amount of charge (whether negative or positive) increases. This 
is also biophysically reasonable since more charge on the impermeable solute 
leads to a greater amount of intracellular countcrions, thus increasing intracel- 
lular osmotic pressure. 

We now turn to the question of stability. Let: 



S = {(c,v) G 



pJV+l 



N 

E 

k=l 



| ^ z k c k + zA/v = 0}. 



(3.18) 



The dynamical system defined by (|3.1[) lives on this set. We must thus modify 
the definition of stability accordingly. A steady state of (|3.1j) is stable if all 
solutions with initial values in S and near the steady state stay close to the 
steady state. A steady state is asymptotically stable if it is stable and if all 
solutions with initial values in S and near the steady state converge to the 
steady state as t — > oo. A steady state is globally asymptotically stable if it is 
stable and if all solutions starting from initial values in S converge to the steady 
state as t — > oo. If we make a change of variables to obtain an iV-dimensional 
dynamical system without the implicit constraint of electroneutrality, the above 
definitions of stability reduce to the usual ones for ODEs. 

We saw in Proposition 12.11 that in the absence of active currents, the free 
energy G defined in (|2.14[) is decreasing. We now construct a free energy like 
quantity that is decreasing in the presence of active currents. Define: 



G = 



N 

•E 

*;=i 



In 



- 1 



% 



A In 



A 



- 1 



(3.19) 



where was defined in (|3.4p . For G, we have the following analogue of Propo- 
sition [57T] 



Lemma 3.2. Let c,v,<j> satisfy system (|3.1[) . We have: 



— - J 

dt ~ ~ ' 

J = (fi + q, L(n + q)) RN + <>w = (l + <lMl + Q) ) m „ +<>w, 



(3.20) 



14 



Mathematics of Pump-Leak Models 



Mori 



where (•, -) RN is the inner product in R . TTie function G is thus a Lyapunov 
function in the sense that it is non-increasing in time. 



Proof. The proof is almost identical to the proof of Proposition (2TTJ The second 
equality in the definition of J comes from (|3.3|) . □ 



If system (|3. 1[) lias a steady state, we may rewrite p.20|) as follows. Let 
(c*, v* , 4>*) be the steady state of (|3.1[) . Note that tt w = at steady state, and 
that qu = —fJ-k where fjj, is the evaluation of the chemical potential at steady 
state. Using this, and the fact that c satisfies the electroneutrality constraint 
(|3.1b|) . we find, after some calculation: 



G(c,v) = G(c,v)-G(c*,v*) 



k=l 



C], I I I \ \ V ) IT 



(3.21) 



The quantity G can be interpreted as being the total free energy of the system 
relative to the steady state. Since G and G differ only by a constant, we may 
replace G with G in (pT20|) : 

dG _ 2 
-jr — — J — ~ \Mi ^Wrn — C^wi 

al . . (3.22) 

p, = Qui, ■ ■ ■ ,mat) t , Jik =ln I -4 ] + Zk((f>- 4>*), 



where we used qk = ^° rewrite J. Thus, if the system has a steady state, 
Lemma 13.21 says that the free energy relative to the steady state is always de- 
creasing at a rate that is controlled by /x, the vector of chemical potential relative 
to the steady state. We shall find both G and G useful depending on context. 

In studying stability, it is sometimes convenient to use a new set of variables 
a = (ai, ■ ■ ■ , ajv) T = vc, v and <fi rather than c, v and (f>. Rewriting (|3.1[) in the 
new variables, a satisfies the differential equations: 

(7a 

- = -L(n + q). (3.23) 



We now state a result on the function G. Note that, although the solutions 
(c,v) of (|3.ip is defined only for (c, v) <G R^ +1 , the function G(c, v) is well- 
defined on R + N x R + where R+ is the set of non- negative real numbers (over line 
of a set will henceforth denote its closure). The same comment applies for G 
viewed as a function of (a, v) and for G. 

Lemma 3.3. 1. Consider the function G defined in (|3.19|) and view this as 
a function of (el, v) G M+ x M+ where a = vc. The Hessian of G is 
positive definite at each point in R^ +1 and G is thus a convex function. 
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2. Suppose condition (|3.5p is satisfied. View G as a function of (c, v) £ 
M + x M + and /ei the unique steady state of (|3.1I) be given by (c*,v*). 
Then, (c*,v*) is the unique minimizer of G restricted to S, where S is 
given in (|3 . 1 8[) . 

Proof. Let Ha denote the (N+l) x(N + l) Hessian matrix of G(a, v). The Hes- 
sian matrix is well-defined for (a, v) £ Kj +1 . For any vector x = (x\,--- ,xn,x v ) £ 



we have: 



X""^ ( 1 v^fe \ 2 ^ 2 
(x,i7 G x) R N + i = > ^=.Tfc + -~x v , (3.24) 

The Hessian matrix Hq is thus positive definite at every point in (a, v) <G 
and G is thus a convex function. 

To prove the second item, we first rephrase the assertion in terms of (a, v). 
We must show that (a*,u*) = (v*c*,v*) is the unique minimizer of G(a,v), 
(a, v) £ R+ x R + when restricted to the hyperplane: 

N 

^2z k a k + A = (3.25) 
fe=l 

where a = (ai, • • • , aAr) T . 

We seek stationary points of G restricted to the hyperplane (|3.25|) . Consider: 

G x = G + \ \J2z k a k + Aj (3.26) 

where A is the Lagrange multiplier. The condition for a stationary point is given 
by: 



0. 



dG x , / a k /v\ 
-g^ = In I — — ) + ZfcA + q k = 0, 

dGx =7] 
dv 



(3.27) 



If we identify A with (f>, the membrane potential, the above condition is nothing 
other than the condition for steady state of system (|3.1[) . Since condition (|3.5[) 
is satisfied, by Proposition 13.11 the above system has a unique solution and 
(a, v,X) = (a*, «*, </>*) where <j>* is the value of (j> at the unique steady state 

of (|3.ip . Since G is a convex function on R + N x R + , its restriction to the 
hyperplane (|3.25[) is also a convex function. Thus, this stationary point is the 
unique minimizer. □ 

We may now state our first stability result. 
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Proposition 3.4. Suppose condition (|3-5[) is satisfied. Then, the unique steady 
state of (|3.1[) is asymptotically stable. Moreover, the decay to the steady state is 
exponential, and the linearized operator around the steady state is diagonalizable 
with real and negative eigenvalues. 

The linearized operator above refers to the linearization when (|3.ip is seen as 
an ODE system on the ./V-dimensional submanifold defined by the electroneu- 
trality constraint p.lbp . We note that asymptotic stability is in fact an im- 
mediate consequence of Lemma 13.31 by a Lyapunov stability argument. Thus, 
if we are only interested in asymptotic stability, there is no need to study the 
linearization. A Lyapunov stability argument will be used to study the global 
behavior of solutions in the proof of Theorem 13.51 

Proof of Proposition \3.4\ It suffices to prove this claim by studying the dynam- 
ics of (|3.1[) in the variables (a, v) instead of (c, v). Since condition (|3.5p is 
satisfied, there is a unique steady state by Proposition 13.11 To study the lin- 
earization around this steady state we must change variables to remove the 
implicit constraint (|3.1b[) (or equivalcntly, (|3.25[0 and obtain an TV-dimensional 
ODE system. 

From (|3.3[) . we have: 



where V a G is the gradient of G(a, v) with respect to a while keeping v fixed. 
Let us examine the matrix L. Take a vector w <E M. . By (|3.3p . we have: 



Since L is symmetric positive definite, the above quantity is non-negative by the 
Cauchy-Schwarzinequality and is equal to if and only if w is a constant multi- 
ple of z. Thus, L is a symmetric positive semi-definite matrix whose eigenspace 
corresponding to the eigenvalue is spanned by z. The restriction of L to the 
orthogonal complement of this eigenspace is thus positive definite. 

Consider a change of coordinates from a to an orthonormal coordinate sys- 
tem b satisfying: 



The JV-th coordinate axis is thus parallel to z. Note that the electroneutrality 
constraint (|3.25[) can be written as 6jv = constant. Let b = £7a, where U is the 
orthogonal coordinate transformation matrix. In the b coordinate system, L 
transforms to L\, = \JLU~ 1 . Given the above properties of L, has the form: 




(3.28) 




(3.29) 



b = (6i, ■ • ■ , b N ) T ', b N = — (z, a) RJV . 



(3.30) 




(3.31) 
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where is a positive definite matrix and Ojv-i is the zero column vector of 
length N — 1. Rewriting (|3.28[) in terms of b, we have: 

db , 
— = -L b V b G (3.32) 

where V&G is the gradient of G (keeping v fixed) seen as a function of b. From 
(|3.31|) . we see that remains constant. Let b = (61, ■ • • , fe^-i). We have: 

^ = -LiV^G (3.33) 

where Vg is the gradient of G with respect to b while keeping b^ and v fixed. 

We have thus reduced the system (|3.1[) to (|3.33[) and to (I3.1cl) which can be 

written as: _ 

dv ,9G , nnt . 
-r = 3.34) 

Letting u = (61, • • • , 6jv-i) v) T , we may write our system as follows: 

1-A.V.g, L.-^ °~-), (3.35) 

where V„ is the gradient of G with respect to u while keeping 6jy fixed. We 
have thus obtained the requisite iV-dimcnsional ODE system in the variables u; 
the electroneutrality constraint 6jv = constant only appears as a parameter of 
the system. 

Let H u be the Hessian of G with respect to u. Given Lemma 13.31 H u 
is symmetric positive definite. Indeed, the quadratic form defined by H u is 
just the restriction of Hq (defined in Lemma 13. 3[) to the subspace of K^" 1 " 1 
orthogonal to (z T ,0) T . The linearized operator of p.35[) around steady state 
is thus given by —L U H* where H* is the evaluation of H u at the steady state. 
Note that L U H* is similar to (iJ*) 1/2 i tl (i?*) 1/2 where (#*) 1/2 is the positive 
square root of H*, which exists thanks to positive dcfinitcness of H*. Since 
(/f*) 1 / 2 i ll (iJ*) 1 / 2 is a symmetric matrix and since L u is symmetric positive 
definite, so is (iJ*) 1 / 2 L„(iJ*) 1 / 2 . Thus, —L U H* is diagonalizablc with real 
negative eigenvalues. The steady state is asymptotically stable and the approach 
to steady state is exponential. □ 

By rewriting (|3.1|l as p.35|) . we see that the system is a gradient flow on the 
hyperplane defined by the electroneutrality constraint where the metric is given 
by L~ . This led us to the conclusion that the linearization is diagonalizable 
with real negative eigenvalues. Note that we made essential use of the symmetry 
of the matrix L, which came from the Onsager reciprocity principle (see (|2.27p ). 
The gradient structure of our system combined with the convexity of G has 
another interesting consequence as we shall see in Lemma [ 
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3.2 Global Behavior 

We now state the main result of this Section. 

Theorem 3.5. Suppose condition (|3.5p is satisfied. Then, (|3.ip has a unique 
steady state and it is globally asymptotically stable. The linearization around 
steady state is diagonalizable with real negative eigenvalues. 

Existence of the unique steady state was proved in Proposition ^. II Asymp- 
totic stability and the property of the linearized operator was proved in Propo- 
sition 13.41 We have thus only to prove that all solutions with initial value in S 
(see p,18[) ) converge to the steady state as t — > oo. Implicit in this assertion 
is that these solutions arc global (defined for all positive time). Once this is 
established, we use the fact that G is a Lyapunov function to obtain the desired 
result. 

To prove that all solutions are global, we must rule out two possibilities. The 
first is that the solution may grow unbounded in finite time. To show that this 
is not possible, we make use of the function G. The second possibility is that 
one or more of the concentrations or the cell volume v may come arbitrarily 
close to in finite time. To show that this cannot happen, we examine the free 
energy dissipation function J defined in (|3.20[) . 

Lemma 3.6. View J defined in p.20|) as a function of (c, v) £ : 

J(c,^J c (c)|(K( M )) 2 , J c (c) = ( 7 + q,£(7 + q)) Rjv - (3.36) 
1. Consider any solution (c(t),v(t)) of system ()3.1[) . We have: 



d 2 A / Pk \" 2A 



p= (pi,--- ,pn) t = £(7 + q)- 



(3.37) 



The function J is thus a Lyapunov function in the sense that it is mono- 
tone non-increasing. 

2. The function J c (c) defined in (|3.36[) tends to +oo as c approaches any 
point on where d- denotes the boundary of a set. 

That J is a Lyapunov function can be seen as follows. View J as a function 
of u introduced in the proof of Proposition 13.41 First, note that: 

¥ = -<V u G,5) = -(VuG,L u V u g) =-J (3.38) 
at \ dt I k jv \ /r n 

where we used (|3.35[) in the second equality. We thus have: 

d J „ / T „ ;x d /„ ~\ \ „ / r „ ~ TT du 



dt 2 { L ^^dA v ^)) P r 2 { L ^^ d t 

= —2 ( L U \7 U G, H U L U S7 U G 



(3.39) 
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where H u is the Hessian matrix of G with respect to u. We used (|3.35[) in the 
last equality. As we saw in the proof of Proposition ^. 41 H u is positive definite. 
Therefore, J is monotone non-increasing. We see that the Lyapunov property 
of J is a general consequence of the fact our system, in suitable variables, is a 
gradient flow of the convex energy function G. 

Proof of Lemma WTb\ We saw above that dJ/dt is non-positive, but we have not 
obtained the right hand expression in (|3.37p . This is most easily done by direct 
calculation. We turn to the second claim. Take any point c b <E and assume 
without loss of generality that the first 1 < I < N components of c b are 0: 



c 



-((),■■■ ,0,c b + i,--- ,c%). (3.40) 



(3.41) 



Decompose the vector 7 + q in the following fashion. 
7 + q= -(7i +7 2 )> 

71 = -{ji + qi,--- ,11 + qi,^_^2) T , 

N-l 

7 2 = -(O,--- ,0,11+1 +qi+i,--- ,iN + qN) T - 

1 

Now, consider a sequence of points c" £ , n = 1, 2, • • • such that c™ — > c b as 
n — > 00. Given that: 

7fe + 9fc=ln^ +q k> k = l,--- ,N, (3.42) 

each of the first I non-zero components of r y 1 goes to +00 as c" — > c b whereas 
7 2 remains bounded as c" — > c b . Now, take an arbitrary vector 

w = (u>i,-- - ,w N ) G R N , |w| = 1, w k > 0, k = l,-- - ,N. (3.43) 

Recall from (13.291) and the subsequent discussion that ( w, Lw) is positive if 

\ /r n 

w is not parallel to z. The vectors w and z are indeed not parallel since z must 
have at least one component that is negative (there is at least one ionic species 
with negative valence). Therefore, 

min (w, Lw) = K w > (3.44) 

w k >0,k=l,- - ,JV,|w| = l \ / R N 

given that the set satisfying Wk > 0, k = 1, ■ ■ • , N, |w| = 1 is compact. There- 
fore, for any vector u whose components are non-negative, we have: 

u,iu\ >K w \u\ 2 . (3.45) 
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Now, let us take the limit of J c (c) as c™ — > c b . If c™ is sufficiently close to 
c b , the first I components of -y 1 as defined in p.41[) are positive. Therefore, we 
have: 



Jc(c n ) > (7i+7 2 ^(7i+7 2 ) 

(3.46) 



> K w | 7l 



Til + (72^72 



where we used (|3.45p and the Cauchy-Schwarz inequality. Since | — > +00 
and 7 2 remains bounded as c™ — > c b , J c (c") tends to +00. □ 

Proof of Theorem \3.5[ Take an arbitrary initial value (c°, v°) £ S where S was 
defined in (|3.18| . We first show that the solution to (|3.1|) starting from (c°,t>°) 
is defined for all t > 0. 

View G of (|3.21| as a function of (c, v). Consider the set: 

Am = {(c,v) e S\G(c,v) < A/}, (3.47) 

where we choose M so that M > G(c°,v°). Given p.22p . the solution stays 
within Am so long as the solution is defined. We first show that the set Am is 
bounded and that it is bounded away from the hyperplane v = 0. 
Any element in Am satisfies: 

It is easily seen that the left hand side of the above is greater than or equal to 
0. Therefore the right hand side must be greater than 0, from which we obtain: 











-n 








v* ) 



(3.49) 

Therefore, v must satisfy 0<w_<w<w+<oo for some constants v + and 
w_. Let M be the supremum of the right hand side of p.48[) over u_ < v < v + . 
This M is clearly finite. Thus, 

N / / / 

Ck 



k=l 



E( c M ln (#) - 1 )+ c l) <M. (3.50) 



We thus see that Ck must be bounded above by a constant c + that depends only 
on M . Therefore, we have: 

< w_ < v < v+, Ck < c + ,k = 1,-- • ,N. (3.51) 

Let (c(t),v(t)) be the solution to (|3 . 1 1) with initial data (c°,v°). Since the 
solution stays within Am, we know that the solution satisfies the bound (|3.5ip . 
We now show that the concentrations Ck(t) are bounded away from 0. Recall 
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from (|3.37p of Lemma 13.61 that the function J is a non-increasing function in 
time. Thus, J(c(t),v(t)) < J(c°,v°) = Mj. Since J c (c) < J(c,v), we have 
J c (c(i)) < Mj. By the second item in Lemma T3.6[ the set: 

{c=(ci,--- ,c N ) eR£|J c (c) <Mj,c k <c+,k = l,---N} (3.52) 

must be bounded away from SM^. Therefore, we have: 

0<c_ <c k (t) <c + ,k = !,-■■ ,N, (3.53) 

where c_ is a constant that depends only on M and Mj. This, together with 
(|3.51|) . implies that the solution (c(t),v(t)) lies in a compact subset JC of S. 
This shows that the solution must be defined for all time. 

We now show that the solution (c(t),v(t)) converges to the steady state 
(c*, v*). Take an arbitrary e > and let B t C R 7V+1 be the open ball of radius 
e centered at (c*,v*). We must show that (c(t),v(t)) £ B t after finite time. 
Observe that we can make 5 > sufficiently small so that As C B e (As is 
defined by replacing M with S in (|3.4T|) ) . This is clear since, by Lemma [3.31 
(c*,v*) is the unique minimizcr of G over S. 

Take S so small that As G B e . If M < S, Am C As C B e . Since the solution 
is contained in Am, there is nothing to prove. Assume M > S. We would like 
to show that the solution is contained in As C B t after finite time. We prove 
this by contradiction. Suppose the solution never enters As- Recall that the 
solution (c(t), v(t)) was contained in a compact set K, C S. The function J(c, v) 
is clearly positive on K\As, since (c*, v*) € As is the only point at which J = 0. 
Since K,\As is a compact set, J > Kj > on K,\As where Kj is a positive 
constant. By (|3.20l) of Lemma 15721 (or equivalcntly, (|3.22l) ). we see that: 

G(c(t), v(t)) < M - Kjt. (3.54) 

This implies that the solution will be in the set As for t > (M — 8)/Kj, a 
contradiction. □ 

Thcorcm l3.5l shows that system (|3.ip has the following remarkable robustness 
property Suppose the pump rates p and the extracellular concentrations c e are 
perturbed within the bounds of condition (|3.5[) . Then, the cell will approach 
the new global steady state. 

The next theorem shows that when condition (|3.5p is not met, the cell volume 
v(t) — > oo as t — > oo. There is thus a dichotomy in the behavior of system (|3.1[) 
depending on whether condition (|3.5[) is satisfied. 

Theorem 3.7. Suppose condition (|3.5[) is not satisfied so that / m i n (q, c°, z) > 0. 
Take any solution (c(t),v(t)) to p. II) . 

1. Suppose / m j n (q, c c ,z) > 0. Then, 

lim v(i) = oo. (3.55) 



22 



Mathematics of Pump-Leak Models 



Mori 



2. Suppose / m i n (q, c e ,z) = 0. Let 4> m i n be as in Provosition \3.1\ Then, 

lim Cfc(t) = c% exp(-q k - z k (j> min ), k = l,---N, 

'7.°° M (3-56) 
lim u(tj = oo. 



Proof. We shall work with the variables c and w = 1/v. In these variables, (|3.ip 
can be written as: 

dc - 

— = w(-L(7 + q) + Ctt w c), (3.57a) 

JV N 

= ^2 z k c k + zAw = (3.57b) 

fc=i fe=i 

w 2 (it w . (3.57c) 



dw 2 



dt 

The solutions are defined on the set: 

jv 

T={(c,w) £ R+ +1 |^z fe c fe + zAw = 0}. (3.58) 

fc=i 

Note that T is just the set S of (|3.18|) written in the (c, w) coordinates. Take 
any initial data (c°, w°) £ T and let (c(t),w(t)) be the solution to (|3.57p starting 
from this point. Showing that v(t) — > oo is equivalent to showing that w(t) — > 0. 
We divide the proof into several steps. 

Step 1: View G defined in (|3.20[) as a function of (c,w). Consider the set: 

JV 

A M = {(c,w) £ M.^ +1 \J2 z kCk + zAw = 0, G(c,w) < M}. (3.59) 

k=l 

This is the same set as (|3.47[) except that we use the function G(c, w) instead 
of G(c,v). We prove that Am is a bounded set. 
For any point in Am, we have: 

J2 (<* (in - 1 + q^j + c%j < Mw - Aw (ln{Aw) - 1) . (3.60) 

Since the left hand side is bounded from below, there is some positive constant 
m, independent of M, such that: 

-m<Mw-Aw(]n(Aw)-l). (3.61) 

Let: 

g ( w ) =-— + AQn(Aw) - 1). (3.62) 
w 
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The function g(w) is a monotone increasing function in w such that g(w) — > — oo 
as w — > and <?(u>) ->ooas«i-> oo. Let w+(M) = g~ 1 (M). Given (|3.61|) . we 
have: 

< w < io+ (M) for any (c, w) G .4m- (3.63) 

Since w is bounded between and w+(M), we see from (13.601) that must also 
be bounded in Am- 

< c fe < c+(M), fc = 1, ■ • - AT for any (c,iy) G ,4m- (3.64) 

S'iep 2: We prove that the solution (c(t),w(t)) is defined for all positive 
time. Choose M = Mo so that Mo > G(cq,wq)- Suppose that the solution 
exists only up to t < To, To < oo. Since G is monotone non- increasing, we have 
(c(t),w(t)) € Amo:^ < To- Since Am is bounded by (|3.63[) and (|3.64|) . there is 
a sequence of times t\ < t 2 ■ ■ ■ To such that x(i n ) = (c(t n ),w(t n )) —> x b = 
(c b , iy b ) S 9^1a/ as n — >• oo. The limit point x b cannot be in R^ +1 since, if so, 
the solution can be continued beyond time To. We also see that Ck, k = 1, ■ ■ ■ ,N 
must stay away from by an argument using Lemma 13.61 similarly to the proof 
of Theorem E31 This implies that c b G and w b = 0. By (|3~Tcj) . we have: 



l(^)=<(|:(«-«+^ 



(3.65) 



Since the right hand side of the above is bounded in Am by (|3.63j) and p.64|) . 
l/w(t) remains finite in finite time. Thus, w(t n ) — > is impossible as t n — ¥ 
To < oo. We have a contradiction. 
Step 3: Let O be the orbit: 

G = {(c(t),w(t))€T,t>0}. (3.66) 

We would like to see whether 

J= inf J(c,w) (3.67) 

(c,ui)GO 

is positive. Recall that J = in T if and only if the point (c, w) is a steady 
state of (|3.1[) . Given our assumption that (|3.5p is not satisfied, by Proposition 
13.11 a steady state does not exist. Therefore, J > in O C T- Thus, if J = 0, 
since O is a bounded set, there is a sequence of points x n G O, n = 1, 2, • • • that 
approaches a point x°° = (c 00 ^ 00 ) G T such that J(x") -> as n -> oo. The 
limit point x°° cannot be in T € R^ +1 since J > there. Since c(t) stays away 
from dRf, c°° <^ 5M+ 7 . Thus, c°° G M+ 7 and = 0. Since the function J is 
continuous up to points (c,w) = (c 00 ,!])^ 00 <G R+, we examine the positivity 
of J on the set: 

K = {x= (c,w) G 5T|c € K^, w = 0}. (3.68) 
On 7?., J can be written as: 

N 

J(c, w = 0) = J c (c) + C(^w(c)) 2 , Tr^(c) = ^ (4 - cfc) . (3.69) 

fc=i 
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We see that J = if and only if J c (c) = and tt^ = 0. It is easily seen that 
J c (c) = in 1Z if and only if c = c* where c* is given by: 

c* = (c l5 - ■ ■ ,c* N ) T , c% = c%cxp(~q k - z k (j) min ), k = l,---N, (3.70) 

where </> m i n is as defined in the statement of Proposition 13.11 Let us evaluate 
7r^ at this point. Substituting (|3.70p into (|3.69p and recalling the definition of 
/min in (O, 

7 r°(c*) = -/ min (q,c e ,z). (3.71) 

Therefore, we have: 

J > if /min > 0. (3.72) 

When / m i n = 0, we have the following. Let B n be the open ball of radius 7/ > 
centered at (c*,0). Then, 

J = inf J(c,w)>0. (3.73) 

' (c,w)eO\B v 

Step 4- We prove our claim when / mm > 0. Given that (c(t), w(t)) £ 0, by 
Lemma 13.21 we have: 

G(c{t), w(t)) < M - Jt. (3.74) 

By p.63p . we have: 

< w(t) < w + (M - Jt). (3.75) 

Note that J > by (|3~?2"j) . Since w + (M) -> as M -> -00, we see from ([3~75|) 
that w(t) — > as t — > 00. This proves (|3.55|) . 

Step 5: In the rest of the proof, we study the / m in = case. As an initial 
step, we prove the following. For any 77 > 0, there is a time t v > such that the 
point (c(t r; ), w^trj)) is in B^. We prove this by contradiction. Suppose otherwise. 
Then, there is an 77 > such that B v n O is empty. 

Wc first show that G is bounded from below in O by a constant M n . Suppose 
otherwise. Then, there is a sequence of points x" = (c n , w n ) £ 0,n = 1,2,--- 
converging to x°° = (c 00 ,?!; 00 ) £ O such that G(c n ,w n ) —> —00. Since G(c,w) 

N 

is a continuous function for c £ R.+ , tj; > 0, the only possibility is that w°° = 0. 
Write G as: 

G(c, 10) = -p(c) + A(\n(Aw) - 1), 
w 

N { { { \ \ \ ( 3 - 76 ) 

p(c)=x;(p fc (in(5)-i+ <ft 



fe=i 



It suffices to show that /o(c°°) > 0. If this is true, we see from (|3.76l) that 
G(c n ,w n ) —> 00, contradicting our assumption that G(c n ,w n ) —> —00. It is 
easily seen by a calculation identical to the proof of Lemma 13.31 that the unique 
minimizer of p(c) under the constraint 



JV 
fe=l 
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is c = c*, at which point p(c*) = 0. Given that (c°°,0) ^ B n , we see that 
c°° ^ c*, and thus p(c°°) > 0. 

By Lemma 13.21 and using (|3.73j) , we have: 

G(c(t),w(t)) <M Q -J v t (3.78) 

so long as (c(t),w(t)) £ B v . Note that l n > by 1(57731 . Thus, if * > (M - 

M v )/J_ v , then G < M^, which contradicts our result that G must be greater 
than M v on O. 

Step 6: We would like to show that there is a positive number rj > such 
that any solution with initial data in B v H T will converge to (c*, 0) as t — > oo. 
If this is true, we can combine this with the result of Step 5 to immediately 
conclude that that all solutions of (|3.57|) converge to (c, w) = (c*, 0) as t — > oo. 
This would prove (j3.56p . 

The vector field defined by the right hand sides of (|3.57a[) and (|3.57cp is 
degenerate at w = 0. Rescaling the vector field by a positive scalar factor does 
not alter the solution orbits, so we shall study the beh avior of an appropri- 
ately rescaled system ( Chiconei . 19991 Benson et al.l . l2010l) . Rescale Q3.57a[) and 



(j3.57c[) by a factor of 1/w. This removes the degeneracy at w = 0: 

dc ~ 

— = -L(7 + q) +(tt w c, (3.79a) 
dr 

^=<7r w . (3.79b) 
dr 

We have taken the time parameter to be r to distinguish the solutions of this 
system with those of (|3.57j) . We consider (|3.79j) on the set: 

N 

T' = {(c,w) elf xR^\J2 z kC k + zAw = 0}. (3.80) 

fe=i 

The difference between T and T 1 is whether or not w is allowed to be equal to 
0. 

If we can show that any solution to ((3.79P with initial data in B v flT converges 
to (c*, 0) as r — > oo, the same is true for p.57|) as t — > oo. This can be seen 
as follows. Let (c(t), w(t)) be the solution to ((3.57P with initial data (c°,w°) £ 
B V DT and (c(r), w(t)) be the solution to (|3.79[) with the same initial conditions. 
Define the function with: 

t= / -—-dr. (3.81) 
Jo w(t) 

This function is well-defined since w(t) > 0. This positivity is a simple con- 
sequence of the backward uniqueness of solutions. The expressions c(£(t)) and 
w(£(t)) satisfy p.57p . and thus, by uniqueness of solutions: 

( C (t), W (t)) = (cm)Mm))- (3.82) 
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By (|3.8ip and the fact that w(t) — > as r — >■ oo, we see that £ — > oo whenever 
t — > oo. Therefore, 

lim (c(i),tu(i)) = Iim (c(r),w(r)). (3.83) 

t— ►OO T — f OO 

S'iep 7: We now show that any solution of (|3.79|) with initial conditions in 
B v n T' converges to (c*, 0) if r} is taken small enough. This will conclude the 
proof. 

Let c(t),w(t) be a solution to (|3.79p . Then, by Lemma T3. 61 we have: 

, N , s 2 

— J(c,w) = -2J2{^=-(^V^) - 2Aw((tt w ) 2 = *(c, w). (3.84) 
T fe=i W c fc / 

Since we have removed one factor of w in (|3.79[) compared with (|3.57|) . one 
factor of w is removed accordingly from the right hand side of (|3.37[) . 

First, we take r > small enough so that J > and * > for£> r n T'\{(c*,0)}. 
By ()3.73p . this is possible for J. Consider If w > 0, the condition ^ = if 
and only if pk — 7r w = 0. This is equivalent to the condition for (|3.57j) to have 
a steady state in T. By Proposition 13. 11 such a point does not exist if / m in = 0. 
For w = 0, the *(c, 0) = if and only if: 

Hi + q) - C< c = 0, (3.85) 

A' 



ZkCk = 0, (3.86) 



fc=i 

where 7r° was given in (|3.69p . It is clear that the point c = c* satisfies the above. 
We show that c = c* is the only point that satisfies both (|3 .85[) and p.86p in a 
hborhood of c = c* in R N . We use the implicit function theorem. Compute 
the Jacobian matrix of the left hand side with respect to c and evaluate this at 
c*: 

B* kl =L kl \+Cc* k . (3.87) 
c i 

Here, is the kl entry of the TV x N Jacobian matrix B* . The rank of B* is 
the same as the rank of B* whose kl entry is given by: 

B* kl =L kl + (c* k ct. (3.88) 

Since L is symmetric positive semidefinite with rank N —1 (sec proof of Propo- 
sition [3T4J) and ( > 0, B* , and thus B* is at least rank N — 1. It is easily 
checked that z c = {z\c*, • • • , zmc* n ) t is an eigenvector of B* with eigenvalue. 
Therefore, all the points that satisfy ()3.85p near c = c* lie on a one-dimensional 
manifold in that is tangent to z c at c = c*. Since (z,z c ) RJV 7^ 0, the only 
common point between this one dimensional manifold and the hypcrplanc (|3.86p 
is c = c* . 
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Define the set: 

T> g = {{c,v) 6 B r DT'\J(c,v) < S}. (3.89) 

Take S = So small enough so that T>s C B r . It is clear by (|3.84[) that any solution 
in T>s stays within this set. Choose r\ small enough so that B V <~)T' C T>s . This 
is clearly possible since J is a non-negative continuous function on B r PI T 1 which 
is only at (c, w) ~ (c* , 0). Take any e < rj. We may choose a Si > so that 
T>s 1 C B e . Given that $ > on T>s \T>s 11 any solution in B v C T>$ will be in 
B e D in finite time. □ 

In the case of / m i n > 0, we do not have a statement on the limiting value of 
c(t) as t — > oo. The limit point (c, w) = (c*,0) may well exist and the limiting 
value c* should satisfy p.85|) . This follows simply by setting the right hand side 
of p.79[) to 0. If there is only one such point, it should be possible to show, 
with the aid of the Lyapunov function J, that this is the single limit point to 
which all solutions converge. This uniqueness, however, is not clear. 

Even though v(t) — ¥ oo as t — ¥ oo regardless of whether / min > or / m i n = 0, 
the rate at which v(t) grows is different. When f m - m > 0, we see, by combining 
(|3.65|) and (|3.75|) and the definition of w+(M) that v(t) grows at most linearly 
with time and faster than any power t a ,a < 1. When / m i n = 0, we expect the 
growth of cell volume to scale like t 1 ' 2 , as can be seen by taking the special case 
p = with initial conditions c = c c . 



3.3 Epithelial Models 

We briefly consider a simple epithelial model. Suppose we have a single layer of 
epithelial cells which separate the serosal and mucosal sides. The concentrations 
of electrolytes in the serosal and mucosal solutions are assumed constant. Let 
these concentrations be denoted c| and c™ respectively. The voltages in the 
mucosal and serosal sides are also fixed at (j) s and </> m . We can write down the 
following model for electrolyte and water balance for an epithelial cell in this 
layer: 

L s fi s - p m - p s , (3.90a) 
Cs<, (3.90b) 

N 

-- Zk4 = 0. (3.90c) 
fc=i 

The definition of the cellular variables c and v are the same as before. In the 
above, L m s are symmetric positive semi-definite matrices and p m ' s are the vec- 
tor of active currents residing on the mucosal and serosal membrane respectively 
which we assume constant. The chemical potentials /x m ' s = (/i™' s , ■ • • , fi^' s ) T 
are given by: 

/C S = ln (^r) +z k (<b-^ s ), (3.91) 



d(vc) 



dt 
dv 

17 = _ Cm7T 



T n m — 



y^zfcCfc 
fe=i 



dt 
zA 
v 



N 
k=l 
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where <j> is the electrostatic potential inside the cell. The osmotic pressure 7r™ ,s 
is given by: 

k=l \k=l ! 

and the hydraulic permeabilities £ m ' 8 are non-negative constants. 

The above problem is in fact mathematically identical to system (|3.1[) . Let: 

C = Cm + Cs, L = L m + L s , (3.93) 
and suppose that £ > and L is symmetric positive definite. Define: 

p = p m +p s -L m /3 m -L m /3 s , 

(m,s\ V * / 

\) + Zk4>m ' & ' 

Then, the triple (c,v,<p) satisfies (|3.90[) if and only if it satisfies (|3.1j) with 
C, L, c c = (cl, ■ ■ ■ , c c N ) T and p prescribed as in (|3.93[) and (|3.94l) . We thus have 
the following result. 

Theorem 3.8. Consider system (|3.90l) . Suppose L m + L s is symmetric positive 
definite and Cm + Cs > 0. Define f mnl as in (|3.5[) in which q = (qi, - ■ ■ , 5jv) = 
L~ 1 p and c%,L and p are prescribed as in (|3.93j) and (|3.94|) . If f m m < 0, t/ie 
conclusions of Theorem \3.5\ hold. If / m i n > ; i/ie conclusions of Theorem \3. 7| 
/io/d. 

Note that this epithelial model also enjoys the robustness property described 
after the end of the proof of Theorem 13.51 We may argue that this is advanta- 
geous for an epithelial cell, which must withstand large changes in extracellular 
ionic concentrations. 



4 Results in the General Case 

In the previous Section, we assumed that the passive transmembrane ionic flux 
jk is linear in fi. In this Section, we establish results that are valid when we 
only assume conditions (|2.24j) . (|2.26|) and (|2.27[) for jk and j w . In particular, 
this will apply to the case when the Goldman equation (|2.3j) is used for jk- We 
also relax the assumption that the pump rates pk be constant. We consider the 
system: 

= -j(<£,M) -ap(^,fi), (4.1a) 

N A N 

= Vv t + z- = Vz fc 4, (4.1b) 
' — ' v * — ' 

fc=i fc=i 

% = -jw(TTw)- (4.1c) 
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The extracellular ionic concentrations c|, k = 1, • • • , N and the amount of im- 
permeable organic solute A are positive. We assume that j, p and j w are C 1 
functions of their arguments. The only difference between this and system 
(|2.11[) is that we replaced p (or p k ) in (|2.f lap with ap (or ap k ) where a is a 
pump strength parameter. We shall find it useful to vary this parameter in the 
statements to follow. 

4.1 Solvability 

We first discuss what we mean by a solution to the initial value problem for (|4.1[) . 
Consider the two constraints (|4.1bl) and (the dimensionless form of) (|2.8[) , which 
we reproduce here for convenience: 

A zA 

Q(c,v)= y]z k c k + — = 0, (4.2) 

k=l 
N 

I{c,(j),a) = ^2z k (j k (c,(t>) +ap k (c,(j))) = (z,j + ap) RN = 0. (4.3) 

k=l 

Define the following set: 

T Q = {y = (c, v,^)elfxi + xl Q(c, v) = I a (c, 4>) = 0}. (4.4) 
We shall often omit the dependence of / and r on a. 

Definition 4.1. Let y° = (c°, v°, (f>°) € V where c° = (c?,--- ,c° N ) T . Let 
c(t) = (ci(t),--- ,CN(t)) T , v(t), t > be C 1 functions and (j>(t), t > be a 
continuous function oft. The function y(t) = (c(t),v(t), 4>(t)) is a solution to 
(|4.1|) with initial values y° if it satisfies (|4.1|) and y(0) = y . The solution may 
or may not be defined for all positive time. 

Since we are solving a differential algebraic system, we must specify ini- 
tial conditions that satisfy the constraints. Note that we require <fi(t) to be a 
continuous function of t. 

Wc have the following on the solvability of (|4.1j) . 

Lemma 4.2. Let yo = (c°,u°,0 o ) G T and suppose dL/d(f> ^ at this point. 
Let B r C M Ar+2 be the open ball of radius r centered at yo. Then, there is a 
K > such that Bk has the following property. 

1. The set Bk DT is an N -dimensional submanifold ofM. N+2 . There is a 
C 1 function $ such that any point y £ Bk H T can be written as y = 
(c,«,$(c)). 

2. Take any point y 1 = (c \v , <j> ) G Bk H T. There is a unique solution 
y(t) = (c(i), v(t), 4>{t)) to (|4.ip with initial values y 1 for short times. For 
short times, <fi(t) = $(c(f)), and thus tfi(t) is a C 1 function. 
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Proof of Lemma \J72] The first item is a straightforward consequence of the im- 
plicit function theorem. For the second item, substitute <j> = $(c) into (|4.1[) . 
Solve this ODE with initial values (c 1 ,?; 1 ) and let c(t) and v(t) be the resulting 
solution. It is clear that (c(i), v(t), $(c(i))) is a solution to (|4.1|) with initial 
values (c l ,v l ,(fr l ) for short times. This solution is the unique solution, since 
4>(t) must be continuous and thus, must remain within Bk for short times. □ 

The same statement clearly holds if we replace Bk with a neighborhood of 
yo. Note that, when z ^ 0, any point in Z^R-nT can be written as (c, V(c), $(c)) 
where V(c) is found by solving Q(c,v) — for v. Thus, when z ^ 0, c serves 
as a local coordinate system on Bk H I\ 

Given the structure conditions on jk we have the following solvability result. 

Proposition 4.3. Suppose jk satisfies (|2.24[) and (|2.27[) . Let 

C r = {y = (c, «, 0) G R N+2 | |c-c c | <r, |0| < r}, (4.5) 

where |-| /or a vector in R N denotes its Euclidean norm. There are positive 
constants K a and K with the following properties. 

1. Take any \a\ < K a . The set Ck H T a is an (unbounded) N -dimensional 
submanifold of M. N+2 such that any y € Ck H T a can be written as y = 
(c, v, <fr(c,a)) where $(c, a) is a C 1 function of c and a. 

2. System (|4.ip ratt initial values y° = (c°,u ,^) ) G Cr- fl r a /ias a unique 
solution y(t) = (c(t), /or s/iori times. For short times, <J)(t) = 
$(c(i),a). 

Proof. To construct the function $(c,a), we use the implicit function theorem 
around c = c c , = 0, a = on /. Note that: 

7(0 = 0, c - c c , a = 0) = (z, j(0 = 0, /x = 0)) RN = (4.6) 

where we used the definition of fi in the first equality and (|2.24l) in the second 
equality. Take the derivative of I with respect to </>: 

d<j> { Z ' d(j> + <9/x Z } RJV + 

In the above, j is viewed as a function of <f> and fi whereas p is viewed as a 
function of <j> and c. We have used the definition of fi to obtain the second term 
in the above. We see that 

dl I d] \ 

— (0 = O,c = c c ,a = O) = (z,-iz) >0. (4.8) 
d(j) \ dn / RN 

where we used ([2725]) . Positivity follows from (f2~27|) . With ({476]) and (gTSJ), we 
can use the implicit function theorem to obtain a C 1 function $ satisfying 

J($(c, a), c, a) = 0, $(c c , 0) = 0, (4.9) 

in a neighborhood of c = c°, a = 0. The rest of the proof is the same as that of 
Lemma 14.21 □ 
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4.2 Existence of Steady States and Asymptotic Stability 

A point y = (c, v, (f>) e T is a steady state of (|4.ip if the right hand side of (|4.1a|) 
and (|4.1c[) is at that point. We have the following result on the existence of 
steady states. This should be seen as a generalization of condition (|1.3I) . 

Proposition 4.4. Suppose j k and j w satisfy (|2.24p . (|2.26[) and (|2.27[) . TTien, 
(|4.1I) /ias a steady state with v > /or a/Z sufficiently small a > so long as the 
following condition is met: 




> 0, 



(4.10) 



0=0,^=0 

where (<9j/<9/x) -1 is t/ie inverse of the Jacobian matrix dj/d/j,. 

The idea behind this result is the following. System (|4.1[) possesses an obvi- 
ous "steady state" when a = 0: c = c c , </) = and u = oo. If a is positive but 
small, we expect this steady state to persist. In order for v to be positive as a 
is perturbed, we need condition (|4.10l) . 

Proof of Proposition \J^4\ Let w = 1/v. Set the right hand side of (|4.1|) equal 
to 0. We have: 

,H) + ap k (<j>,fi) = 0, 



3k( 



1, 



,N, 



N 

E 

fc=i 



ZkCk 



zAw 



Jw (^*w ) 



0, 
0. 



(4.11) 



View the above as an equation for c, <fi and w. Note that c = c c , (f> = 0, w = is 
a solution to the above system if a = 0. To apply the implicit function theorem, 
we show that the Jacobian matrix with respect to c, (f>, w is invertible at this 
point. This is equivalent to showing that the only solution to the following 
linear equation for c = (ci, • • ■ , c/v), 4> and w is the trivial one. 



N 



\ - djk_ 



djv 











0=0, /i=0 \ I 





0, k = 1, 



,iV, 




zylw = 0, 



AQ = 0, 



(4.12) 
(4.13) 
(4.14) 



where we used (|2.24j) (and its consequence (|2.25j) ') to obtain (|4.12j) . Equation 
(|4. 12[) together with condition (|2.27|) and (|4.14|) together with condition (|2.26p 
gives: 

N 

0, ^S + Au; = (4.15) 

fe=i 



c fc - 
-r + Zfc<; 
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where the first equation is holds for all k. Using (|4. 1 5[) to eliminate w and c k 
from (|4.13p . we have: 



JV 

E 

fe=i 



N 



(z k - z)z k cl<j) = zl c k$ = 



(4.16) 



k=l 



where we used (|4.1b|l in the first equality. Since c k is positive and at least one of 
z k 7^ 0, we see that = 0. From (|4. 15[) . we see that c k — for all k and w = 
since A > 0. We can thus invoke the implicit function theorem to conclude that 
we have a solution c(a), <fi(a) and w(a) to (|4.11j) when a is close to 0. To ensure 
that w (or cquivalcntly, v) is positive for small a > 0, we compute: 



dw 
da 



Wc»,ffL 



(4.17) 



i)=0,fi=O 



Since w = at a = 0, condition (|4.10p will ensure that the v is positive for a 
small and positive. □ 



We may also compute d<j)/da: 



da 



N 

E 
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(zc c - z c ), 



a=0 \ k= i / \ 

z c = {zic\, ■ ■■ , z N c e N ). 



dj_ 



i=0,^t=O 



(4.18) 



Given (|4.10p . this shows that the sign of 4> is the same as the sign of z if \z\ is 
large enough. In physiological situations, z is large and negative, and thus we 
will have a negative membrane potential. 
Condition (|4.10j) applied to (fTTT]) yields: 



3fNa^ 



2[K" 1 



>0, 



5Na 



9K 



(4.19) 



thus reproducing condition (|1.3[) . It is interesting that we do not have any 
restriction on z for this to be true. Given (|4.17p . the cell volume will be small 
if (|4.10p is large. Expression (|4. 19[) is indeed large: for a typical cell, [Na + ] c ^> 
[K+] e and g Na < g K - 

We now turn to the question of stability of steady states. 



Definition 4.5. Suppose y* = (c*,v*,(f>*) € T is a steady state of (|4.ip . Let 

B r be the open ball of radius r centered at y* . The steady state y* is stable if 
the following is true. For any small enough e > ; there exists a S > with the 
following property. Choose anyy G BsCiT . Then, the solution(s) y(t) to system 
(|4.1I) with initial value y° is defined for all positive time and satisfies y(t) € B e . 
The steady state is asymptotically stable if it is stable and all solutions with 
initial data y° £ T sufficiently close to y* approach y* as t — > oo. 
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Since any solution lies on T, we may replace B e with B t n T. The only 
difference between the usual definition of stability and the one above is that the 
initial data must lie on T. If dl /dip ^ at y* = (c*, v* , <fi*) and z ^ 0, then, by 
the remark after the proof of Lemma 14.21 system (|4.1[) can locally be written 
as an ODE for c only. The above definition of (asymptotic) stability is then 
equivalent to the (asymptotic) stability of c* for this ODE system. 

Let c* = (c*, • ■ ■ , c* N ), v*, cf)* be a steady state of (|4.1[) . Define the following 
quantities: 




(4.20) 



If j w satisfies f|2 . 23[) . .7 w (7r w ) = if and only if ir w = and thus tt^ = 0. In this 
case, 7r w = 7r w . 

Let G be free energy with respect to the steady state defined in p.21|) . We 
have the following analogue of Proposition 12. II or Lemma I3~2l 

Lemma 4.6. Suppose c* = (c*,--- , c* N ) T , (j>*, v* is a steady state of (|4.ip . 
Then, we have: 

dG N 

fe=i 

where jk — jk — ik > Pk = Pk — Pk an d 3k > Pk are the passive and active fluxes 
evaluated at the steady state. 

Proof. Rewrite (|4.1a[) as follows: 

9<yV Q^ = + a Pk) + tit + apt) = -^2% (j k + apt^j , (4.22) 

fc=i 

where we used j| + ap* k = 0. Note that, j^, the water flux at steady state, 
is equal to 0. Thus j w = j w — = j w . The rest of the proof is the same as 
Proposition 12. II □ 

If we apply the above lemma to system (|3.1j) . this is nothing other than 
(|3.22|) . The next Lemma gives us a sufficient condition for asymptotic stability 
in terms of G. 

Lemma 4.7. Let y* — (c*, v*,<f>*), c* = (c*,--- ,c* N ) T be a steady state of 
(|4.ip . Suppose there is neighborhood U C K. Ar+2 of y* such that U H T is an 
N -dimensional submanifold in which any point y G IA H T can be written as 
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y = (c, w, $(c)) where <£> is a C 1 function of c. Suppose any solution y(t) = 
(c(i), v(t), 4>(t)) inlA (or equivalently, inUHT) satisfies: 

^<-A^(|£| 2 + |£ w | 2 ) (4.23) 

for some positive constant K*. Then the steady state is asymptotically stable 
and the approach to the steady state is exponential in time. 

Proof. As in the proof of Proposition 13. 4[ we will find it convenient to use the 
variables a = (aj,--- ,ajy) T — vc,v and <f> rather than c, v and <fi. We shall 
continue to use the symbols U,T to denote the corresponding sets in the new 
coordinates. View G as a function of a and v. Note first that: 

G(a>*) =0, 

= (7fe -7fc)l c =c* =0, 

a=a.,„=„. (4-24) 
= 7r w| c=c * /u=1 ,» = 0, 

a— a* ,v—v* 

where a* = (aj, • • • , a^) T = v*c* . By Lemma T3.3[ G is a globally convex func- 
tion on (a, v) <G R^ +1 . The point (a, v) = (a*,v*) is thus the global minimizer 
of G. The positive definiteness of the Hessian matrix of G implies that there is 
a neighborhood N C R N+1 of (a*, v*) where 

K G l (|a-a*| 2 + \v-v*\ 2 ^) < G(a,v) < K G (ja - a*| 2 + |v - v*\ 2 ^j (4.25) 

for some positive constant Kq. 

Now, consider R^" 1 " 2 with the coordinates (a, v, <f). Define Q ~ Y^,k=i z k{ a k~ 
Of,). This Q is the same as the Q in (|4.2[) except that it is written in terms 
of a. We claim that, in the vicinity of (a* ,v* ,<fi*) in 1 N+2 , the set of vari- 
ables (£t, 7? w , Q) defines a coordinate system. It is easily seen that the variables 
(c,i>,</>) defines a coordinate system. We thus consider the coordinate change 
from (/i,7f w ,(5) to (c, v, <fi). The Jacobian matrix between these two sets of 
variables at (c,v,(j>) = (c*,v*,(j>*) is non-singular. This computation is almost 
the same as the computation in the proof of Proposition 14.41 so we omit the 
details. The claim follows by the implicit function theorem. There is therefore a 
neighborhood V C R*^ 2 of (a* ,v* ,4>*) in which the following inequality holds: 

K, (|a-a*| 2 + \v-v*\ 2 + \<f>- 0*| 2 ) < \fi\ 2 + l^w| 2 + \Q? (4.26) 

where is a positive constant. Any solution to (|4.1I) satisfies Q = (see (14.21) ). 
Thus, we have: 

K^(\a-a*\ 2 + \v-v*\ 2 ) < \fi\ 2 + \ir w \ 2 (4.27) 
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for any solution in V. 

Choose a neighborhood M C R N+1 of (a*,v*) such that (a, v, <&(a/v)) g 
U n V for all (a, n) e M. Consider the following set: 

G = {(a, v) e M Ar+1 | G(a, w) < M G , M G > 0} (4.28) 

Since G(a.,v) is a convex function such that G(a*,v*) — 0, Q is an open neigh- 
borhood of (a*,v*), and can be made arbitrarily small by making Mq small. 
Take Mq so small that Q C Af n A4 . Consider the following open neighborhood 
of (a*, v*,^): 



W = {y = (a, v, 4>) e R w+2 | (a, »)e5,yeWn V}. (4.29) 

Any solution in W, by definition, belongs to W PI T. For any such solution, we 
have: 



_<„^(|£|2 + |£ w | 2 ) 

< (|a - a*| 2 + |„ - v*\ 2 ) < -^^G 



(4.30) 



where we used (|4.23l) in the first inequality, (|4.2T[) in the second inequality and 
(|4.25[) in the third inequality. Solving the above differential inequality, we have: 

G < Maex^i-K^K^t/Ka), t > 0. (4.31) 

We thus see that W fl T is a positively invariant set, and thus, all solutions 
starting from W n T are defined for all time. By (|4.25[) . (a.(t),v(t)) approaches 
(a*,v*) exponentially in time. Since (WnT) C (UnT), (f> = $(a(t)/w(t)). Since 
$ is a C 1 function, <j)(t) also approaches qb* exponentially in time. □ 

We are now ready to state the main result of this Section. 

Theorem 4.8. Suppose jk and j w satisfy (|2.24[) . (|2.27[) and (|2.26|) and jk,Pk 
and c c k satisfy (|4. 10|) . For all sufficiently small a > 0, the steady states found 
in Proposition are asymptotically stable. The approach to steady state is 
exponential in time. 

In Proposition 13.41 we used the symmetry condition of (|2.27|) to show that 
the eigenvalues of the linearized matrix around steady state are all real. Here, 
we cannot prove such a statement. In fact, the proof to follow goes through 
even if we assume f|2 . 28[) instead of (|2.27p . 

We also point out that, unlike Proposition 13.41 or Theorem l3.5[ we can only 
draw conclusions when the pump rate is small (a is small). One may wonder 
whether it may be possible to generalize Theorem 14.81 to the case when the 
pump rate is not necessarily small. For this, one would clearly need a condition 
stronger than (|2.26[) or (|2.27[) . One natural idea is to require that jk and j w 
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satisfy (|2.26p and (|2.27p not only at tt w = and fi = but at any arbitrary 
value of 7r w and [i: 



-(ttw) > for all 7r w , 
0, jLt) is a positive definite matrix for all and 



(4.32) 



This stronger condition is indeed satisfied by the Goldman equation (|2.3p (and 
trivially by (|2.2[1 ). Let us assume the pump rates apk are constant. A natural 
conjecture may be that if j w and jk satisfy (|4.32j) . any steady state of (|4.1[) is 
stable (whether or not the pump rate is small). This statement is true provided 
that jk is only a function of it and not a function of 0. In the case of (|2.2j) or 
(|3.ip . jk is indeed only a function of fi. Unfortunately, (|2.3p is a function of both 
/x and 0. The danger when jk depends on independently of fi is that djk/dcj) 
may adversely affect the stability properties imparted by condition (|4.32l) . 

When the pump strength is small (a small), djk/d(j) is small provided a is 
small thanks to condition (|2.25|) . This is one of the key observations that we 
will use in the proof to follow. 



Proof of Theorem [7T5[ Let y* = (c*,v*, <f>*) be the steady state found in Propo- 
sition 331 As a > is made small, y* £ Ck H L dchncd in Proposition 14.31 We 
shall henceforth assume that y* 6 Ck H T. 

By Proposition ^. 31 Lemma |4~B1 and Lemma |4~71 it is sufficient to show that, 
for sufficiently small a, there is a neighborhood U C R JV+2 of (c*,v*,<fi*) such 
that any y = (c, v, <fi) £ U n L satisfies the following inequality: 

N 

K * (iai 2 + i^wi 2 ) < ^2 ^ k {p k + a P k ) + ^ wjw = J ( 4 - 33 ) 
fe=i 

for some > 0. The neighborhood U may depend on a. 

Recall from the proof of Lemma 14.71 that (J1,tt w ,Q) defines a coordinate 
system in the vicinity Af of (c* , v* , (j)*) . Note that the point (c*,v*,4>*) is the 
origin in the coordinate system (/2,7r w ,Q). Let: 

V r = {(c,v,cf>) GAT | \£\ 2 + 9l + Q 2 <r 2 ,r > 0} (4.34) 

and take r small enough so that T> r C Ck- Define the set: 

r Q = {y = (c, v, 0) e R+ x R + x M | Q(c, w) = 0}. (4.35) 

Since T C Tq, it is clearly sufficient if we can find a small enough r > such 
that (|4.33p holds for any point in V r n Lq . 

We can write = 0—0* as a function of /2, 7r w and Q in 2? r . In particular, 
on T> r H Lq, is a function of p, and 7r w only. Call this function = (p(p,Tr w ). 
The function <^ satisfies: 



fc - z)c]J (exp (jufe - z fe ^) - 1) - Z7fw- (4.36) 
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This equation is obtained by expressing c k and v in terms of fl k , 7? w and <fi and 
substituting this into Q = 0. 

Now, take any point w = (/I 1 , 7?^, 0) =€ 2? r nrg, /2 1 = (juj, ■ ■ ■ , ^jy) T where 
we have expressed the point w using the (/2,7r w ,(3) coordinate system. Let us 
compute the right hand side of (|4.33[) at this point. 

N 

fc=i 

(4.37) 

where we took j k , p k as functions of 4>, M and j w as a function of 7? w . For j k , 
we have: 

j k ((p(ti 1 ,^l f ),fi i ) = J -^j-jk (^p(s'fi 1 ,STf w ),sli 1 ^ ds 

= % n IA{%Wt + w) {s ^^A ds (4 - 38) 



f 1 fdjk dtp ! \ 



Performing a similar calculation for and j w and substituting this back into 
(I4.37|) . we obtain the following expression. 

= {ip,Pip) MN+1 , if} = (fi 1 ,^), 

f 1 , , (4-39) 

P= / (L + B + C)(sfi\ S nl)ds : 
Jo 

where L, B and C are (AT + 1) x (N + 1) matrix-valued functions defined on 
£Vnr<g, given as follows. Let L^, £?m and Cki be the fcZ entries of these matrices. 

'dh/fffii ifl<k,l<N, 
Lki = { dj w /dn w if k = I = N + 1, (4.40) 
otherwise, 

\dj k /d$){dv/djli) if I <k,l<N, 

{(&3 k /d$)(d<p/87? w ) i£l<k<N,l = N + l t (4.41) 
otherwise, 

'a(dp k /&jii + (8$ k /d$)(d<p/8fii)) if 1 < ft, I < N, 
Cm = { a(dp k /d0)(d<p/dn w ) if 1 < ft < N, I = N + 1, (4.42) 

otherwise. 

To show that (|4.33[) is valid in T> r n Tq, it is sufficient to show that L + B + C 
is positive definite in T> r n Tq in the sense that: 

(x, {L + B + C)x) RN+1 > |x| 2 (4.43) 
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for any x £ M. N+1 with a constant > that does not depend on the point 
in T> r n Tq. Since L, B and C are continuous functions on V r n Tq and we may 
take r as small as we want, all we have to show is that (|4.43|) holds at the origin, 
(/x, 7r w , Q) = (0, 0, 0), or equivalently, at the steady state. 

Let L*,B* and C* be the evaluation of the three matrices at steady state. 
Since the steady state is a function of a, L*, B* and C* are functions of a. We 
will show that (|4.43|) holds for sufficiently small a > 0. 

First, let us examine the behavior of /x*, ir^, </>* as a function of a. By 
Proposition 14.41 c*, <p* are C 1 functions of a that approach c e , respectively as 
a — > 0. Therefore, /x* is a C 1 function of a that approaches as a — > 0. It is 
clear that 7r^, = for any a. 

The fcZ entry of the matrix L* is given by: 



T* - 
^ki — 



(djk/dm)^,^ itl<k,l<N, 
(dj w /d-K w )\ 



if fc = 2 = TV 
otherwise. 



(4.44) 



as a — > 0, given (|2.27|) and (|2.26p . there is a constant 



Since /x* — s- and <^ 
Id, > such that: 

(x,L*x) RN+1 >K L \x\ 2 (4.45) 

for sufficiently small a > 0. 

To examine B* and C*, let us first compute dip/dj2i and dtp/dTT w . This can 
be computed by taking the partial derivatives of (|4.36|) : 



dip 
dfii 

dip 



(Zl - z)Cj 



£=O,7f w =0 



^=O,7f w =0 




(4.46) 



Since c* k 



c| and 1/v* — > as a — >■ 0, we see that both of the above are 
bounded (and in fact has a definite limit) as a — > 0. 

Let us examine B*. For 1 < k, I < N, the fcZ entry of the matrix B* is given 

by: 

l 

djk 



B* kl = (zt 




(4.47) 



0=0*, M=M* 



Recall that djk/dcj) = at /x = from (|2.25p . Since </>* — > and /x* -> as 
a -> 0, by (|2.24j> . we see that -B^ -> as q -> 0. Likewise, -B^ — > as a -> 
when 1 < k < N and I = N + 1. 

Now, consider C^, the kl entries of the matrix C* . For 1 < k, I < N, we 
have: 

-l 



I ^Pk . I \ * 




dpk 



(4.48) 



,/J,=H' 
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Given that /x* > 0, (f>* — > as a — > 0, the quantity inside the outer-most 
parentheses remains bounded as a — > 0. Thus, — > as a — > 0. The same 
conclusion holds for C^, 1 < fc < AT, Z = AT + 1. 

Since L* satisfies (|4.45[) from sufficiently small a and J3* and C* both tend 
to the zero matrix as a -> 0, L* + B* + C* satisfies (|4.43[) with A'* = K L /2 for 
small enough a. □ 



5 Discussion 

In this paper, we presented what the author believes is the first analytical re- 
sult on the stability of steady states of pump-leak models. In Section [3J we 
studied the case in which the flux functions are linear in the chemical poten- 
tial. In the proof of Proposition 13. 4[ we saw that the system can be seen as a 
gradient flow of a convex function. This is nothing other than the relaxation 
law postulated in linear non-equilibrium thermodynamics, dX./dt = —LVxG, 
where X is the vector of extensive variables, L is the matrix of transport co- 
efiicients and V xG is the gradient o f the free energy G with resp ect to X 



( Katzir-Katchalskv and Currant Il965t iKielstrup and Bedeauxl . 120081 ). There 



are two interesting points here. The first point is that this gradient flow is 
restricted to flow on a submanifold on which electroneutrality holds. The elec- 
trostatic potential, as we discussed in the proof of Lemma |3~31 can then be seen 
as a Lagrange multiplier of this gradient flow. The second point is that we can 
find a suitable modification of the free energy (G or G) so that our system is 
a gradient flow even in the presence of active currents. We proved the follow- 
ing results. If condition (|3.5[) is satisfied, there is a unique steady state that is 
globally asymptotically stable. If not, the cell volume tends to infinity as time 
t — > oo. The system is thus robust to external perturbations in the following 
sense. Suppose the cell is subject to a change in extracellular concentration or 
pump rate that stays within the bounds of condition (|3.5[) . The cell will even- 
tually approach the new global steady state. We saw that the same conclusions 
hold for a simple epithelial model, since it could be mapped to the single cell 
problem. 

In Section [H we proved that steady states for pump- leak models are stable 
so long as the steady state is not too far away from thermodynamic equilib- 
rium. The "stable equilibrium state" is the "state" at which all intracellular 
ionic concentrations Ck are equal to the extracellular concentration c| and the 
cell volume v is infinite. If the pumps work in the "right direction" (in the 
sense of Proposition I4.4[) v can be made finite even with a small pump rate. A 
biophysical interpretation of Theorem l4.8l is that if the pump rate is sufficiently 
small, the new steady state of finite volume is still close enough to thermody- 
namic equilibrium so that the steady state inherits the stability properties of the 
equilibrium state. It is interesting that stability of thermodynamic equilibrium, 
which may be considered the "dead" state, confers stability to the "live" state. 

The results of Section [4j though applicable to general pump- leak models, 
only asserts the existence of at least one asymptotically stable steady state for 
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small pump rates. To draw analytical conclusions at large pump rates without 
the linearity assumption of Section [3l it is likely that one would need to look at 
special characteristics of specific pump-leak models. Our current study may thus 
be seen as complementing computational investigations of s tability, in which one 



is not constrained to small pump rates ( Weinstein . 19971 ). We also point out 



that, for general pump-leak models, we cannot rule out the possibility of multiple 
stead y states or of other non-trivial asymptotic behavior. Indeed, ([Weinstein . 
19921) " reports an instance in which there are two stable steady states in a non- 
electrolyte model of epithelial cell volume control. 

Many epithelial models include effects not included in model (12 . 1 IP or p. 901) . 
Of particular impor t ance i s the incorporation of acid-base reactions (| Wein stein. 
19831 IStrieter et al.l . 1990). It is usually assumed that the acid-base reactions 



are sufficiently fast so that these reactions are in equilibrium. This gives rise to 
additional algebr aic constra i nts, increasing the co-dimension of the d ifferential 
algebraic system (IWeinsteinl . I2002L l2004t IWeinstein and Sontad.120091) . It would 
be interesting to see whether the analysis of this paper can be extended to this 
case. A starting point for an analytical study of such models will likely be a 
free energy identity. A potential complication is that the algebraic constraints 
of acid-base reactions are not linear in the concentrations. This may pose diffi- 
culties in extending the global results of Section [3J The author hopes to report 
on such an extension in future work. 

Stability of steady states is ju st a starting point in the study of homeostatic 
contr ol in epithelial systems. In ( Weinsteinl . |2002[ 12004 ; Weinstein and Sontad . 
20091 ). the authors go beyond stability to study the optimal control of home- 



ostatic parameters by minimizing a quadratic cost function along a relaxation 
trajectory. We hope that our current study will lead to new insights into such 
problems. 

Free energy dissipation identiti es similar t o (12.131) are present in ma ny models 
of soft-condensed matter physics ( Doi and Edwards . 1988 ; Doi . 20091 ). and can 
be used as a gu i ding principle in formu lating models in dissipative systems 



( Eisenberg et all l20ld iMori et all l201lh . The present work owes much of its 



inspiration to this body of work. We believe that there is much to be gained 
by a systematic application of these ideas to the study of physiological systems. 
We hope that this paper will be a starting point in this direction. 
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